********************************************************************************
* 'Revisiting the effect of GnRH analogue treatment on bone mineral density in 
*  young adolescents with gender dysphoria'
* Dr Michael Biggs, Department of Sociology, University of Oxford
* Journal of Pediatric Endocrinology and Metabolism, 2021
********************************************************************************

capture: program drop sumofsquares
program define sumofsquares, rclass

syntax , p(real) mean(real) n(int)

local t = invttail(`n'-1, `p'/2) * cond(`mean'<0, -1, 1)
di "t = " `t'
local sd = (`mean' * sqrt(`n') ) / `t'
di "sd = " `sd'
local ss = `n' * ( (`sd')^2 * ((`n'-1)/`n') + (`mean')^2 )
di "Sum of squares = `ss'"
return scalar sumofsquares = `ss'

end

*****************************

capture: program drop aggregatesamples
program define aggregatesamples, rclass

syntax , p_A(real) mean_A(real) n_A(integer) p_B(real) mean_B(real) n_B(integer)

di "Sample A"
sumofsquares, p(`p_A') mean(`mean_A') n(`n_A')
local SS_A = r(sumofsquares)
di _n "Sample B"
sumofsquares, p(`p_B') mean(`mean_B') n(`n_B')
local SS_B = r(sumofsquares)
di _n "TOTAL"
local n_T = `n_A' + `n_B'
di "n = `n_T'"
local mean_T = (  (`mean_A'*`n_A') + (`mean_B'*`n_B')  ) / `n_T'
di "mean = `mean_T'"

local sd_T = sqrt(   (  ( (`SS_A'+`SS_B') / `n_T' ) - (`mean_T')^2 )  *  (`n_T'/(`n_T'-1))  )
di "s.d. = `sd_T'"

local t = `mean_T' / (`sd_T' / sqrt(`n_T'))
di "t = `t'"

local p = ttail(`n_T'-1, abs(`t') ) * 2
di "p = `p'"
return scalar p = `p'

end

********************************************************************************

capture: program drop analyze_density
program analyze_density
syntax name, figurenumber(integer) bonelabel(string) bonelabelshort(string) //, varlabel

local variable "`namelist'"
di "`variable'"

local varlabel = strproper(subinstr("`namelist'", "haz", "", .)) // ("`namelist'", 1, strpos("`namelist'", "haz")))
di "`varlabel'"

ttest `variable'zb=`variable'z12 if t==0
ttest `variable'z12=`variable'z24 if t==0
ttest `variable'z24=`variable'z36 if t==0

gen `variable'z = `variable'zb if t==0
replace `variable'z = `variable'z12 if t==1
replace `variable'z = `variable'z24 if t==2
replace `variable'z = `variable'z36 if t==3

summ `variable'z
local max = max(r(max), abs(r(min))) +.1
local min = -1 * max(r(max), abs(r(min))) -.1
di `min'
di `max'
 
local p05 = invnormal(.05)
di `p05'
local p01 = invnormal(.01)
local p001 = invnormal(.001)
di `p001'
local p0001 = invnormal(.0001)
local z2 = round(normal(-2)*100,.01)
di `z2'
local z3 = round(normal(-3)*100,.01) 
di `z3'

gen `variable'z_under2 = cond(`variable'z==., ., `variable'z < -2)
tab `variable'z_under2 if t==2
sort `variable'z
list `variable'z if t==2 & `variable'z < -2, clean noobs

sort `variable'zb t
twoway ///
	(connected `variable'z months if id==1) ///
	(connected `variable'z months if id==2) ///
	(connected `variable'z months if id==3) ///
	(connected `variable'z months if id==4) ///
	(connected `variable'z months if id==5) ///
	(connected `variable'z months if id==6) ///
	(connected `variable'z months if id==7) ///
	(connected `variable'z months if id==8) ///
	(connected `variable'z months if id==9) ///
	(connected `variable'z months if id==10) ///
	(connected `variable'z months if id==11) ///
	(connected `variable'z months if id==12) ///
	(connected `variable'z months if id==13) ///
	(connected `variable'z months if id==14) ///
	(connected `variable'z months if id==15) ///
	(connected `variable'z months if id==16) ///
	(connected `variable'z months if id==17) ///
	(connected `variable'z months if id==18) ///
	(connected `variable'z months if id==19) ///
	(connected `variable'z months if id==20) ///
	(connected `variable'z months if id==21) ///
	(connected `variable'z months if id==22) ///
	(connected `variable'z months if id==23) ///
	(connected `variable'z months if id==24) ///
	(connected `variable'z months if id==25) ///
	(connected `variable'z months if id==26) ///
	(connected `variable'z months if id==27) ///
	(connected `variable'z months if id==28) ///
	(connected `variable'z months if id==29) ///
	(connected `variable'z months if id==30) ///
	(connected `variable'z months if id==31) ///
	(connected `variable'z months if id==32) ///
	(connected `variable'z months if id==33) ///
	(connected `variable'z months if id==34) ///
	(connected `variable'z months if id==35) ///
	(connected `variable'z months if id==36) ///
	(connected `variable'z months if id==37) ///
	(connected `variable'z months if id==38) ///
	(connected `variable'z months if id==39) ///
	(connected `variable'z months if id==40) ///
	(connected `variable'z months if id==41) ///
	(connected `variable'z months if id==42) ///
	(connected `variable'z months if id==43) ///
	(connected `variable'z months if id==44) ///
	, ///
	yscale(range(`min' `max')) ///
	ylabel(-3(1)3, labsize(vlarge) angle(0) nogrid)  ///
	ymtick(-2 -3, notick grid glcolor(black) glwidth(medthick) glpattern(dash)) ///
	ytitle("{it:z}-score", size(vlarge)) ///
	xlabel(0 12 24 36, labsize(vlarge)) ///
	xtitle("Months on GnRHa", size(vlarge)) ///
	plotregion(margin(l=0 t=0 b=0 r=2)) ///
	legend(off) graphregion(fcolor(white) lcolor(white) margin(zero)) name(bone, replace)

graph twoway function y=normalden(x,0,1), range(`min' `max') lw(medthick) lcolor(blue) horizontal color(blue)  /// recast(area)
	ytitle("") ///
	ylabel("") ///
	ytick(-2 -3, glcolor(black) glwidth(medthick) glpattern(dash)) ///
	yscale(range(`min' `max')) ///
	xlabel(0 "", labsize(vlarge) labcolor(none)) ///
	xtitle("Normal", size(vlarge)) ///
	plotregion(margin(zero)) ///
	fxsize(20) ///
	legend(off) graphregion(fcolor(white) lcolor(white) margin()) ///	
	name(normal, replace)

count if `variable'z!=. & t==0
local n = r(N)
graph combine bone normal, imargin(0 0 0 0)  ///
	title("{bf:Fig. `figurenumber'.} `varlabel' `bonelabel' over time ({it:n} = `n')", ///
		size(large) color(black) margin(b=3)) ///
	graphregion(fcolor(white) lcolor(white) )
graph export "Figures/Fig_`variable'_time.jpg", as(jpg) replace


count if `variable'z!=. & t==2
local n = r(N)
summ `variable'z if t==2
local mean = r(mean)
local sd = r(sd)
graph twoway ///
	(histogram `variable'z if t==2, width(1) start(-4) freq fcolor(gs8) lcolor(gs8)) ///
	(function y=`n'*normalden(x,0,1), range(-4 4) lw(medthick) lcolor(blue) color(blue)) ///
	, ///
  title("{bf:`varlabel' `bonelabelshort'}", color(black) margin(b=3) span) ///
  xtitle("") ytitle("") ///
  xlabel(-4(2)4, labsize(large)) ///
  xmtick(-4(1)4) ///
  yscale(on) ///
  xscale(range(-4 4) noextend) ///
  ylabel(0 5 10,labsize(large) angle(0)) ///
  ytitle("") ///
  xscale(range(-4 4)) ///  lw(medthick) 
  plotregion(margin(zero)) ///
  legend(off) ///
  graphregion(fcolor(white) lcolor(white) margin(r+3)) ///
  name(`variable', replace)

end

********************************************************************************
********************************************************************************
********************************************************************************

cd "/Users/michael/Documents/Gender/Health/Tavistock/Bone density/"

use "../EarlyInterventionData/854413_GD_PublicDataset_v1.0.dta", clear
* Data from https://reshare.ukdataservice.ac.uk/854413/

save scratch, replace

********************************************************************************
* Table 1

* Note where the table gives p = .000 (to 3 decimal places) then p must be smaller than .0005; use that as upper bound 

local n_M = 10
local hip_mean_Mb =    .45 
local hip_mean_M24 =  -.600
local spine_mean_Mb =   .130
local spine_mean_M24 = -.890
local spinehaz_mean_Mb =   .486
local spinehaz_mean_M24 = -.279
local hip_M_p =      .002
local spine_M_p =    .0005 // upper bound
local spinehaz_M_p = .0005 // upper bound

local n_F = 21
local hip_mean_Fb =  -1.075
local hip_mean_F24 = -1.779
local spine_mean_Fb =  -0.715
local spine_mean_F24 = -2.000
local spinehaz_mean_Fb =  -0.361 
local spinehaz_mean_F24 = -0.913
local hip_F_p =      .001
local spine_F_p =    .0005  // upper bound
local spinehaz_F_p = .001 

gen Xrow = "Measure" in 1
replace Xrow = "Subsample" in 2
replace Xrow = "Mean Z-score at b" in 3
replace Xrow = "Mean Z-score at 24" in 4
replace Xrow = "Change in Z-score" in 5
replace Xrow = "p-value" in 6
replace Xrow = "N" in 7

foreach measure in hip spine spinehaz {
	gen X`measure'_Josephetal = cond("`measure'"=="spinehaz", "spine BMAD", "`measure' BMD") ///
			if Xrow=="Measure"
	gen X`measure'_2011to15 = ""	
	replace X`measure'_Josephetal = "Joseph et al." if Xrow=="Subsample"
	replace X`measure'_2011to15 = "2011-15" if Xrow=="Subsample"
	foreach t in b 24 {	
		local `measure'_mean`t' = (``measure'_mean_M`t'' * `n_M' + ``measure'_mean_F`t'' * `n_F') / (`n_M' + `n_F')
		replace X`measure'_Josephetal = string(round(``measure'_mean`t'', .01)) ///
				if Xrow == "Mean Z-score at `t'"
		quietly summ `measure'z`t' if `measure'z24!=.0
		local `measure'_2011to15_`t' = r(mean)
		replace X`measure'_2011to15 = string( round(r(mean), .01) ) ///
				if Xrow == "Mean Z-score at `t'"
		}
	foreach sex in M F {
		local `measure'_mean_`sex' = ``measure'_mean_`sex'24' - ``measure'_mean_`sex'b'
		}
	replace X`measure'_Josephetal = string( round(``measure'_mean24' - ``measure'_meanb', .01) ) ///
			if Xrow == "Change in Z-score"
	replace X`measure'_2011to15 = string( ( round(``measure'_2011to15_24' - ``measure'_2011to15_b', .01) ) ) ///
			if Xrow == "Change in Z-score"
	quietly aggregatesamples, mean_A(``measure'_mean_M') n_A(`n_M') p_A(``measure'_M_p') ///
							  mean_B(``measure'_mean_F') n_B(`n_F') p_B(``measure'_F_p')
	replace X`measure'_Josephetal = string( round(r(p), .001) ) ///
			if Xrow == "p-value"
	quietly ttest `measure'zb=`measure'z24
	replace X`measure'_2011to15 = string( round(r(p), .001) ) ///
			if Xrow == "p-value"	
	replace X`measure'_Josephetal = string( `n_M' + `n_F' ) ///
			if Xrow == "N"
	replace X`measure'_2011to15 = string( r(N_1) ) ///
			if Xrow == "N"
	}
export excel X* using table1.xlsx if Xrow!="", replace
drop X*

/*
* Test of formula
local p_A = 0.0042201
local mean_A = -6.750
local n_A = 4
local p_B = 0.0571910
local mean_B = 2.6666667
local n_B = 3
aggregatesamples, mean_A(`mean_A') n_A(`n_A') p_A(0.002) ///
	mean_B(`mean_B') n_B(`n_B') p_B(0.001) 
*/

********************************************************************************
* Figures

expand 4, gen(dup)
by id, sort: gen t=_n-1
list id t
gen months = t*12+runiform(-1, 1)

analyze_density hip, figurenumber(1) bonelabel("bone mineral density") bonelabelshort("BMD")
analyze_density spine, figurenumber(2) bonelabel("bone mineral density") bonelabelshort("BMD")
analyze_density spinehaz, figurenumber(3) bonelabel("bone mineral apparent density") bonelabelshort("BMAD")



foreach variable in hip spine {
	count if `variable'z!=. & t==2
	local n_`variable' = r(N)
	}

graph combine hip spine spinehaz, ///
	imargin(0 1 3 0) xcommon ycommon hole(2)  /// 
	title("{bf:Fig. 1.} Bone density after 24 months of puberty suppression", ///
		size(large) color(black) margin(b=3)) ///
	l1title("Frequency", size(medsmall)) ///
	b1title("Z-score", size(medsmall)) ///
	graphregion(fcolor(white) lcolor(white) ) ///
	caption("N = `n_spine' for spine, `n_hip' for hip. BMD, bone mineral density; BMAD, bone mineral apparent density.")	
graph export "Figures/Fig_combined.jpg", as(jpg) replace

list spinez if t==2 & spinez < -3, clean noobs
di "Proportion where Z below -3 is " (round(normal(-3)*100,.01)) "%"


